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Abstract 


The second version of the computer program DRENAS (subsurface 
agricultural drainage) is presented, which includes four analytical 


solutions and a finite element solution Galerkin type of the differential 
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equation of agricultural drainage in its one-dimensional form and can be 
run under Windows systems of 32-bit or 64-bit. The analytical and 
numerical models included in version 1.0 of the program were migrated 
from Visual Basic 2005 to Visual Basic 2017 and modified to expand and 
improve their calculation, data recording and results visualization 
capabilities. The numerical model can now solve the non-linear 
Boussinesq equation with variable coefficients including the vertical 
recharge term and allows use at the boundary of the drains, Dirichlet 
condition, linear radiation or fractal radiation. The results provided by 
the module of the numerical solution of the program were validated 


using an analytical solution and a drainage test. 


Keywords: Boussinesq equation, variable storage coefficient, fractal 


radiation condition, soil hydraulic properties. 


Resumen 


En este trabajo se presenta la segunda versión del programa de 
cómputo DRENAS (drenaje agrícola subterráneo), el cual incluye tres 
soluciones analíticas y una solución de elemento finito tipo Galerkin de 
la ecuación diferencial del drenaje agrícola en su forma unidimensional; 
puede ejecutarse bajo sistemas operativos Windows de 32 bits y 64 bits. 
Los modelos analíticos y numéricos incluidos en la versión 1.0 del 
programa se  reprogramaron en Visual Basic 2017, realizando 
modificaciones para mejorar sus capacidades de cálculo, registro de 
datos y visualización de resultados. El modelo numérico se amplió para 
resolver la ecuación de Boussinesq no lineal con coeficientes variables, 


incluyendo el término de recarga vertical, y para manejar las opciones 
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de usar en la frontera de los drenes condiciones tipo Dirichlet, radiación 
lineal o radiación fractal. Los resultados proporcionados por el módulo 
de la solución numérica del programa fueron validados 
satisfactoriamente haciendo uso de una solución analítica y un 


experimento de drenaje. 


Palabras clave: ecuación de  Boussinesq, coeficiente de 
almacenamiento variable, condición de radiación fractal, propiedades 


hidráulicas del suelo. 
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Introduction 


Agricultural land that presents problems of salinity and/or shallow water 
table can be recovered or controlled by artificial drainage, which in the 
case of the farm can be of the subsurface type. The variables to be 
determined in the design of a subsurface drainage system are the depth 
of the drains, separation and diameter. In calculating these variables, 
the characteristics of the soil (texture, structure, physical and hydraulic 


properties) and the condition of the water flow in the soil must be 
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considered. 


It is possible to carry out studies of water transfer processes in 
agricultural drainage systems with the Boussinesq equation of the 
unconfined aquifers, which although it simplifies mass and energy 
transfers in the vadose zone of the soil, provides general descriptions of 
the water flow in the saturated zone (Ritzema, 2006). For detailed 
modeling of water dynamics with the Boussinesq equation, in the case of 
unconfined aquifers, the dependence of the storage coefficient on the 
hydraulic head should be considered and the boundary condition in the 
drains that best represents the water flow drained should be used. On 
the one hand Fuentes, Zavala € Saucedo (2009) formally establish the 
relationship between the soil-water retention curve and the storage 
coefficient in unconfined aquifers, and on the other hand Zavala, 
Fuentes 8 Saucedo (2007) demonstrate that the transfer of water from 
the soil into the drains must be described with a non-linear radiation 


condition. 


Zavala, Saucedo 8 Fuentes (2014) develop the first computational 
version to apply the Boussinesq equation with variable storage 
coefficient subject in the drains to fractal radiation conditions which they 
call DRENAS (subsurface agricultural drainage). However, the software 
was programmed in Visual Basic 2005, so it operates only on 32-bit 
Windows Systems and its graphic and database options depend on 
Microsoft Office 2003. 


The aim of this paper was to develop in Visual Basic 2017, the 
second version of the DRENAS computer program to expand its 


calculation capabilities, improve its graphical interface and eliminate its 
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dependence on external programs for data management and graphics 
development. In the second version of the software, new calculation 
options were incorporated and programmed to allow simulation 
alternatives such as the handling of the recharge or discharge term of 
the Boussinesq equation as a time-dependent function, representing the 
initial condition of the hydraulic head as a variable function in the space, 
possibility of using mechanistic restrictions for the shape parameters of 
the Van Genuchten model for the soil-water retention curve (Van 
Genuchten, 1980) when used to define storage coefficient as a function 
of hydraulic head and the option of use on the drains Dirichlet boundary 
conditions. With the purpose of debugging the computer program, its 
validation was carried out considering an analytical solution for 
subsurface drainage under steady state and data from an experimental 


laboratory drainage test. 
Materials and methods 


Base equations 
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In the computer program, the complete or simplified form of the one- 
dimensional Boussinesq equation of agricultural drainage is solved. In 


general it can be written as: 


mE = Arm] + Ro (1) 


where H = H(x, t)= hydraulic head measured from an impervious 
layer or a reference level (L), and is a function of the horizontal 
coordinate x and time t; T(H) = aquifer transmissivity (L?T*), in the 
case of unconfined aquifers it is T(H) = KsH; Ks = saturated hydraulic 
conductivity (LT?); u(H) = storage coefficient that can be a function of 
the hydraulic head H; Ry(t) = is the volume of recharge or discharge in 


the unit of time per unit area of the aquifer (LT?). 


The recharge or discharge Ry(t) is a difficult variable to determine 
since it depends on the flow conditions existing in the soil surface (rain, 
evaporation, water over soil), the water flow in the soil and the humidity 
regime in the vadose zone, the stratigraphy of the unsaturated zone of 
the soil, macroporosity, cracks, etc. In this new version of the numerical 
simulation model the term R,(t) was incorporated representing it by a 


time-dependent polynomial: 


Ry (t > 0) = a, t? + b,t? +c,t + d, (2) 
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where a,, b,, C, and d, = coefficients to be calculated considering 
experimental measurements or estimates of the recharge time evolution 


for a specific event. 


The relationship between the storage coefficient and the soil-water 


retention curve established by Fuentes et al. (2009) was considered: 
u(H) = 0, — 0(H — H;) (3) 


where 8; = saturated volumetric water content (L%L”*); and Hs = 
elevation of the soil surface (L). It is possible to obtain an explicit 
analytical representation of the storage coefficient from equation (3) if 
the soil-water retention curve is known. The classic van Genuchten 
relationship (van Genuchten, 1980) was selected, because it is widely 
accepted in field and laboratory studies because of its descriptive 


flexibility : 
0(p) = 0, + (05 — 0,11 + (Y/Pa)"]"" (4) 


where 0, = residual volumetric water content (1913); y <0 = scale 
parameter of the soil-water pressure potential vw, (L); m and n = 
dimensionless shape parameters. This new version of the computer 
program included mechanistic relations between m and n that satisfy 
the infiltration theory, the classic Burdine restriction (Burdine, 1953) " m 
= 1- 2/n" as well as the fractal relationships of Fuentes, Brambila, 


Vauclin, Parlange € Haverkamp (2001) that derive in their study of the 
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hydraulic conductivity of unsaturated soils. Fractal constraints are the 
so-called neutral pore " m = (1 - 4s/n)/s ", geometric pore " m = (1 -— 
2s/n)/s " and big pore " m = (1 - 4s/n)/2s "; where s = quotient 
dimension of the fractal object, as the soil is considered in that study, 
"s" being the ratio between the fractal dimension of the object " Ds" 
(Mandelbrot, 1982; Falconer, 2014) and the Euclidean space "E", s = Df 
/E. 

By entering equation (4) into equation (3) the following 


relationship was obtained: 
u(H) = (0s — 0,1 — [1+((H = Hs) /payiI”) (5) 


The fractal radiation condition of Zavala et al. (2007) was retained 


in this second version of the computer program: 
a = YKinl(Ha — D,/P1% (6) 


where qa = drain flow; y = dimensionless conductance coefficient; 
Kin = WKsKg) = soil-drain interface conductivity, K; = saturated hydraulic 


conductivity of the soil, Ky = drain wall conductivity; P = depth of 


drains; Ha = hydraulic head in the drain position; Do. = aquifer 
thickness; s = (1/2)(Si:+S4q) = quotient dimension of the soil-drain 
interface; si= quotient dimension of the soil; and sy = quotient 


dimension of the drain wall. The relationship between areal porosity 
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(farear) and volumetric porosity (dp) is obtained according to the 
probabilistic idea as Harea, = dp *= y, it is necessary that volumetric 


porosity is: 


(1-P+p7=1 (7) 


and areal porosity: 


1 
1 = 
(l — Marcar): + eee =1 (8) 


The conductivity of the drain wall is obtained with a formula based 
on Poiseuille's law: Ky = (1/21(9/W)area ARuo)?, where g = gravity 
acceleration; v = water kinematic viscosity; Ryp = hydraulic radius of 
the drain. The option to directly manage the linear radiation condition (s 
= 0.5 in equation 6) was introduced in version 2.0 of DRENAS without 


using the s;, and sy routines. 


Simplified numerical modeling of agricultural drainage can be done 
by imposing Dirichlet type boundary conditions (hydraulic head in the 
drain position). In this new version of the computer program, this 
boundary condition was incorporated as an alternative to the radiation 


condition (6); the next function that depends on the time was selected: 
Hart > 0) = aAagt + bat? + Ca + dut (9) 
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where Xdren = horizontal coordinate where the drain is located; y 
dar Da, Ca and dy = coefficients that must be calculated from 
measurements taken of the evolution of the hydraulic head on the drain. 
Equation (9) has the flexibility to describe extreme behaviors of the 


hydraulic head above the drain (rise and fall water level). 


Finally, to perform the numerical simulation with equation (1), the 
initial state of the hydraulic head in the system must be defined, in this 
work the option of using up to a third degree polynomial that is a 


function of the horizontal coordinate x is included: 
H(x,t= 0) = apx? + byx? + cpx + dy (10) 


where ap, bp, Cp and do, = coefficients to be determined from 
measurements of the hydraulic head in the system at the initial moment 


of calculation or simulation. 


In the numerical solution of the system (1)-(10) the Galerkin type 
finite element method was used for spatial discretization, a finite 
difference scheme for temporal discretization, Picard's iterative method 
for linearization of the resulting system and a Preconditioned conjugate 
gradient method for the solution of the system of algebraic equations 
(Zavala et al., 2014; Zienkiewicz, Taylor, £ Zhu, 2013; Noor € Peters, 
1987). The system of equations resulting from discretization was 


programmed in Visual Basic 2017. 
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It is known that simplified calculations of separation between 
drains and groundwater abatement can be performed by applying 
analytical solutions obtained for reduced forms of the Boussinesq 
equation (equation 1). For example, for the steady flow regime there is 
the Hooghoudt formula (Hooghoudt, 1940) and in unsteady flow the 
relations of Glover-Dumm (Dumm, 1954) and Fuentes et al. (1997), 
which were already included in the original version of the computer 


program. 
Graphical interface development 


In the second version of the DRENAS computer program, all the 
information capture modules of the original version were reprogrammed 
in Visual Basic 2017, independent internal databases of Microsoft Access 
were developed, the processing of the graphics was disconnected from 
Microsoft Excel, and the executable file was generated to install DRENAS 
2.0 on any computer that has a 64-bit Windows Operating System and 


even the option for 32-bit Windows Systems was also generated. 


The DRENAS 2.0 program has two calculation modules, one called 
analytical solutions and another numerical model (Figure 1). The 
analytical solutions module has four sections, two to perform 


simulations for steady water flow conditions using the Hooghoudt 
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solution (Hooghoudt, 1940) and two for unsteady flow conditions, one to 
perform calculations with the Glover-Dumm solution (Dumm, 1954) and 
the other with the solution of Fuentes et al. (1997). In these four 
sections it is possible to calculate the separation between drains or the 
drainage module. Regarding the original version of the DRENAS model, 
in version 2.0 these sections were improved by incorporating the 
following alternatives: a) Storage of simulations performed in an 
internal database; b) Loading of previous examples; c) Generation of a 
simulation report and alternative to export it to a pdf file; d) The results 
graph can also be exported as an image. An example of a calculation 


section of the analytical solutions module is presented in Figure 2. 


VO 1 1 10 0 0 E 1 O 1 0 1 0 0 1 0 O IIA 


Figure 1. General screen of the Figure 2. A calculation section 
DRENAS 2.0. of the analytical solutions 


module. 


The numerical solution module was programmed in a Tab Index 
form that contains six sections: four are developed for the capture of 
data related to the drainage system, drain characteristics, soil properties 


and data necessary for numerical simulation such as those related time 
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control, spatial discretization, boundary conditions and graphical display 
of results; The fifth tab serves to store the sample data in an internal 
database and the sixth section contains the button to execute the 


numerical simulation. 


The new capabilities of the numerical simulation module that were 


programmed were established as follows: 


a) First, from the UNSODA 2.0 database (Nemes, Schaap, Leij, 8 
Wósten, 2001), soil-water retention curves were selected (208 
soils), were adjusted with Equation (4) considering the 
mechanistic restrictions between m y n described above (Zavala, 
Saucedo, € Fuentes, 2018) and the results obtained were 
included in an internal database of the computer program (Figure 
3). If the option to use data from this program base is selected, 
when choosing the type of soil required, the sections of physical, 
hydraulic properties and storage coefficient are automatically 
filled with their corresponding parameters. The improvement to 
highlight in this new version is that the relationships between m 
and n of Burdine models, neutral pore, geometric pore and big 


pore can be managed (Figure 4). 
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Figure 3. Numerical module: expanded UNSODA base. 
LD 
Figure 4. Numerical module: soil-water properties. 
b) In the numerical solution of the Boussinesq equation, the 


recharge term (R,(t)) was incorporated as a function of time 
(equation 2), the necessary programming was carried out to be 
able to use the Dirichlet type boundary condition (equation 9) and 
the initial condition described by equation (10) was also 
programmed, and in Figure 5 the new section developed to 
include these three conditions is presented. 


Figure 5. Numerical module: Limit conditions and recharge. 
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Cc) A new section was developed to record the simulation data, load 
data to compare simulations as well as to select the type of graph 
to be displayed in real time that does not depend on external 
programs (Figure 6). 


Figure 6. Numerical module: simulation graphs. 


Results and discussion 


The results of the computer program were compared against external 
results (analytical solution and a laboratory drainage test) to detect and 


correct programming errors and check the consistency of the coded 
solutions. 
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Validation 1 


An analytical solution of the Hooghoudt type presented in Fragoza et al. 
(2003) that considers linear radiation (s = 0.5 in equation 6) in the 
drains to describe the evolution of the free surface in an subsurface 
drainage system with constant recharge (R¿) and steady-state water 


flow: 
HG) = HH (BR KIMAL=2) (11) 


where H¿ = hydraulic head to center between drains. If defined 


h(x,t) = H(x,t) - Do and consider the linear radiation condition: 


ha = [y (+ y)?D2 +8(Q + y)h¿(h¿ + 2D,) — (4 + y)Do]/[2(2 + y)] (12) 


where ha = hydraulic head above the drain (x=0 and x=L); Do = 
aquifer thickness; and y = conductance coefficient of the linear radiation 
condition used in the drains. Analyzing equation (2) it is deduced for this 


case that a,=b,=c,=0 and d,=Ro. 
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The data of the drainage system are those reported by Fragoza et 
al. (2003) and correspond to a field experiment conducted in the 
irrigation district 076 Valle del Carrizo, Sinaloa. The drainage system 
evaluated has a separation between drains L = 50 m, drain depth P = 
1.5 m, depth of the impervious layer Do, = 3.5 m, elevations of the 
impervious layer and the soil surface H;=0.0 m and Hs=5.0 m. The soil is 
a Clay with a saturation hydraulic conductivity of Ks = 0.557 m/d and a 
storage coefficient y = 0.1087. It is assumed y = 1.5 and h¿ = 0.5 m, 
and with equation (12) allows to obtain hy = 0.365 m and the recharge 
value obtained with the condition for steady-state flow that the drain 
discharge is equal to the recharge is Ra = 0.000944 m/d. With these 
values, equation (11) was applied to obtain the spatial distribution of 


the hydraulic head along the separation between drains. 


The described scenario was reproduced with the computer 
program, modeling the unsteady process of water flow in the system 
until reaching the steady state, which corresponds to that described with 
the analytical solution. Initial condition of constant hydraulic head was 
assumed H(x, t=0) = 4.50 m (ap = bp= Cp= O and dy, = 4.50 m in the 
equation 10) and to have the same conditions of the analytical solution, 
the linear radiation boundary condition was used in the program, s =0.5 
in the equation (6). The hypothesis of equality between the conductivity 
of the soil and the conductivity of the drain wall was considered K; = Kg, 
which implies Kin = Ks = Kg. The conductance coefficient value reported 
by Fragoza et al. (2003), to adapt it to the form of radiation that is 
handled in DRENAS 2.0, getting y =0.045. Spatial discretization was 
performed with the number of nodes N = 1000 (Ax = 0.05 m) and in the 
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discretization of time was used Atin; = 1.157E-05 d, Atmax = 6.94E-04 d 
and Atmin = 1.157E-06 d. The simulation time was 720 days. 


This comparison is evidence that the information capture modules 
of the numerical solution of the program and in particular the calculation 
module that handles the term of the vertical recharge, initial condition 
and boundary condition, operate correctly and are free of errors of 
programming; it is also observed that the simulated series are stable 


and free of numerical oscillations (Figure 7). 


a) First week of drainage b) One year of drainage 


Figure 7. Simulation of an agricultural drainage scenario (constant 


storage coefficient and constant vertical recharge). 


Validation 2 


Tecnología y ciencias del agua, ISSN 2007-2422, 11(2), 262-290. DOI: 10.24850/j-tyca-2020-02-07 


ES 2020, Instituto Mexicano de Tecnología del Agua 
Tecnología y 


CienciaszAgua Open Access bajo la licencia CC BY-NC-SA 4.0 


(https://creativecommons.org/licenses/by-nc-sa/4.0/) 


Zavala, Fuentes € Saucedo (2004) conducted an experiment in a 
drainage module in which two drains of 30 cm in length (/p) and 5 cm in 
diameter (Dp) were installed, with N, = 233 circular perforations of 
diameter d. = 0.158 cm evenly distributed on the surface of each tube. 
The dimensions of each drain: L = 100 cm, P = 120 cm and Do = 25 cm. 
The module was filled with sandy textured soil that was saturated by 
applying a constant water head on the soil surface until the air was 
removed. With the drains covered, excess water was removed on the 
surface and covered to eliminate evaporation. Finally, the drains were 
uncovered and the volume of drained water measured for 10 days (240 
hours); it should be noted that the initial condition corresponds to 
H(x,0) = P + Do (d, = 145 cm in the equation 10 and the rest of its 


coefficients are null) and the recharge is null (R,, = O in the equation 2). 


Soil porosity was 0d, = 0.5396 cm*/cm? and with this value the 
equation (7) was solved obtaining si; = 0.7026. The saturation hydraulic 
conductivity value is K¿ = 18.3 cm/h. The aquifer storage coefficient is 
described with equation (5) for the parameters 0; = 0;,, 0, = 0, Wa = - 
41.8 cm and n = 2/(1-m) = 3.19 (Burdine, 1953). 


The porous areal of the drain ¡iS Harea y = Ao/Aa = 0.0097 cm*/cm? 
and when solving the relationship (8) Sy = 0.5688. With the hydraulic 
radius of the circular holes in the drain wall Ryp = 0.0395 cm, water 
kinematic viscosity v = 36 cm*/h, Poiseuille's law is applied to calculate 
the conductivity of the drain wall K¿ = 2,670.8 cm/h. The conductivity 
values in the soil-drain interface and the relative dimension are Kin = 
221.08 cm/h and s=0.6357. The coefficient y of fractal radiation 
reported by Zavala et al. (2004) was reinterpreted to adapt it to the 
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form of radiation that is handled in DRENAS 2.0 (equation 6), obtaining 
y = 0.0749. 


The variable storage coefficient and the initial condition of total 
saturation of the porous medium cause the simulation to be executed by 
starting small time steps at the beginning of the simulation, which 
minimizes problems of convergence and numerical stability. The time 
steps used are Atini = 2.77E-04 h, Atin; = 2.77E-05 h and Atmax = 0.2 h. 


The number of nodes to discretize the space is 200. 


In Figure 8 (a and b), the simulation data capture process and the 
evolution of the drained depth are presented. In Figure 8b it can be 
seen that the computer program successfully reproduces the 
experimental data of the drained depth, having a stable and oscillation- 
free numerical solution, which is an indicator that the variable storage 
coeficient module works correctly. The time and space steps used in the 
simulation are adequate since the calculated total drained depth (23.915 


cm) differs from the measured depth (23.92 cm) in only 0.03%. 


a. y e FrmGrafica =D ER 


a) Information capture: variable b) Depths drained: experimental 
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storage coefficient. and simulated. 


Figure 8. Simulation of an agricultural drainage scenario (variable 
storage coefficient, null recharge and fractal radiation condition in the 


drains). 


Conclusions 


The new version of the DRENAS computer program (version 2.0) is a 
useful tool to describe the hydraulic operation of subsurface drains, to 
design a drainage system and solve inverse problems of hydrodynamic 
characterization of soils from drainage events. This version of the 
computational tool incorporates the following functions: a) internal 
bases for storing the data of the examples analyzed; b) expanded 
UNSODA (soil database) with parameters for hydraulic conductivity and 
soil-water retention models of the neutral pore, geometric mean pore 
and big pore; b) option to represent the aquifer recharge or discharge 
as a function of time; d) option to represent the initial condition of the 
hydraulic head as a function of space and to use Dirichlet boundary 
conditions that vary over time; e) generates simulation reports and 
allows export of results graphs. The graphical display of the modeling 
avoids dependence on external programs and runs on 32-bit and 64-bit 


Windows Operating Systems. 
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